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Computational methods are developed to determine the proton 
energy degradation and flux attenuation as a function of penetration 
depth in various materials. The primary purpose of this work is to 
estimate the energy deposition or tissue dose rate at a given depth 
or on the surface of a shielded object. Detailed analysis of the methods 
used and their accuracy are a prime part of this study. Numerous 
comparisons are made with results of other workers in this field. 

/jo T6 o^ 


TABLE OF CONTENTS 


Page 

I. INTRODUCTION 1 

II. ENERGY SPECTRA OF PRIMARY PROTONS 1 

III. MULTILAYER SHIELDS !S 

IV. NONELASTIC PROTON COLLISIONS AND 

SECONDARIES 22 

V. PROTON DOSE RATE CALCULATIONS 34 

VI. ISOTROPIC INCIDENT PROTON FLUX ON SLABS 53 

VII. CONCLUSIONS 71 


i i i 



NASA -GEORGE C. MARSHALL SPACE FLIGHT CENTER 


TECHNICAL MEMORANDUM X-53063 


THE CALCULATION OF PROTON 
PENETRATION AND DOSE RATES 


by 


Martin O. Burrell 


RESEARCH AND DEVELOPMENT OPERATIONS 
RESEARCH PROJECTS LABORATORY 





I. INTRODUCTION 


There have been several calculational methods developed to 
determine the proton energy degradation and flux attenuation as a function 
of penetration depth in various materials, the ultimate purpose being 
to estimate the energy deposition or dose rate at a given depth or on the 
surface of a shielded target such as a man. The methods range from 
fairly simple approximations to complex and tedious numerical methods. 
However, most of the methods are essentially the same in that they 
assume the so-called "straight -ahead model. " In this mode 1 , the 
assumption is made that energetic protons lose energy by ionization 
losses associated with the removal of bound electrons in the shield 
materials , 1 with no subsequent change in particle direction. Elastic 
scattering is assumed to be strongly in the forward direction with a 
negligible energy loss and hence is ignored as a slowing -down mechanism. 
However, in most of these models, an attenuation correction is made 
for non -elastic collisions that completely remove the primary proton. 

The degree of sophistication in the non-elastic collision calculation is 
a function usually of the shield thickness and the subsequent treatment 
of the secondary particles liberated. 

The methods introduced by the writer are in the same category 
as those discussed above. It is hoped that the innovations presented 
will help in obtaining reliable solutions in a simpler manner than is now 
available . 


II. ENERGY SPECTRA OF PRIMARY PROTONS 

It seems to follow that regardless of the methods or models used, 
the slowing-down energy loss of the primary protons is assumed to be 
dependent only on the ionization loss from bound electrons 2 which is 
given by various modifications of the Bethe -Bloch formula for stopping 
power; 


1 Hydrogen shields should probably be excepted. 

2 

An additional discussion on this, point may be found in Section IV of this 
report, "Non-elastic Proton Collisions." 
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where E is the kinetic energy of the proton, Z is the atomic number, 
V = ]3 C is the proton velocity, m is the electron mass, N is the 
number of atoms of the material per cm 3 , I is the average ionization 
potential of the material, and C is a correction term for electron- 
shell binding. 


A quantity of greater utility in many of the computational scheme s 
is the range of a proton which is given by 

R ' E > * f il?) ■ (2 > 

o 

The dimensions of stopping power, S(E), are usually (Mev-cm 2 /gm) 
and therefore for the range the dimensions are (gm/cm 2 ) from Eq. (Z) . 

In order to see how the above quantities enter into the calculation of 
proton penetration, the following development is presented. Figure la 
illustrates the parameters of the problem, where E denotes the inci- 
dent energy of a proton and E the energy at depth X . 

Now if certain liberties are granted it can be seen that the proton 
energy in going from E to E' 1 ' might be represented by an analytical 
relationship such as 

E = g(E*) , (3) 

where, obviously, E' 1 ' is a function of X. Hence, the proton differential 
energy flux at depth X > 0 may be related to the flux at depth X=0 by 
a simple change of variable technique denoted by 


^ vt* 

<t> x (E ) dE"' 




dg(E*) 
dE * 


dE' 


( 4 ) 


0£ course, the practicality of representing the flux at depth X , 
as shown in Eq. (4), depends on the ability to find a usable relationship 
between the energy E and E"\ However, the ability to write Eq. (3) 
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FIGURE lb. VARIATION OF SPECTRUM SHAPE AS PROTONS PENETRATE A SHIELD 


in a nice mathematical expression does not follow from direct applica- 
tion of the Bethe -Bloch formula. In order to arrive at a practical 
solution to the problem, one can resort to the following exercise in 
functional manipulation. 

The proton range is assumed to be represented by an empirical 
curve fit, or even as a tabulated set of numbers, in the case of a pure 
numerical approach. Thus, if 


R = F Z (E) [gm/cm 2 ] (5) 

is used to denote the range of a proton of energy E incident on a ma- 
terial denoted by the subscript Z , then at the depth X (gm/cm 2 ) in 
material "Z", the energy of the proton is reduced by an amount AE 
associated with an equivalent reduction in range given by AR = X. 
Thus we can write 


R - X = F z (E - AE) . (6) 

Now E - AE = E'*' , the energy of the proton at depth X , and since 
R = F Z (E) we write 

F Z (E) = X + F Z (E*), and 
E = g(E*) = F [ F z (E*) + X] . f 7) 

Thus, Eq. (7) provides the relationship required by Eq. (3). However, 
there are some obvious restrictions to the functional form which the 
approximation of R(E) can assume. For this reason, use is often 
made of the numerical approaches to finding the proton differential energy 
flux at a depth X . However, it should go without saying that the num- 
ber of functional forms which are amenable to the manipulations indi- 
cated in Eq. (7) are, mathematically speaking, without limits. The 
most popular attempt to arrive at a simple solution to the proton pene- 
tration problem is that given by assuming the range of a proton in a 
material n Z M can be represented simply by 

R = aE r , (8) 

where the coefficient "a" is dependent on the material, and the power 
"r" only slightly dependent on the "Z n number. 1 In fact, a value of 


l - This type of approximation dates back to 1947. R.R. Wilson, Phys. 

Rev., 71, 385L, Chap. 22, Sec. 3 (1947). 


4 



r = 1.78 will suffice for Z = 6 to 30. This choice of range formula 

is usually considered valid from about 1 0 to 250 Mev with a maximum 
error of ± 5% in approximating the various numerical integrations for 
range based on the Bethe -Bloch formula for stopping power. As an 
illustration of the techniques that can be used to arrive at a simple 
formula for primary proton penetration the following is presented: 

Assume that the incident proton energy spectrum is given by 




HE 


-q 


< 


protons 
crn^ -Mev 


E< E; 


(9) 


and that for the slab thickness and energy spread the range is sufficiently 
well approximated by Eq. (8); then, from Eq, (7), we write 
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from which it is readily seen that 


r X / r 

= (E - - ) if E 
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From Eq. (11), it follows that if the slab thickness is exactly X = aE* , 
the incident proton of energy E Q will just reach zero energy at depth X. 
Next we find 


dg(E*) 

dE^ 


E 


r - 1 




r -1 
r “ 


( 12 ) 


Substituting the appropriate results of Eqs. ( 9) —(12) into Eq. (4) we 
obtain 
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</> X (E*) = 
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E* r + — 
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where Eq. (11) must be satisfied for the limits. Figure lb depicts the 
general appearance of the transformations between Eqs. (9) and (13). 


Equation (13) gives the proton differential energy spectrum at 
depth X for the incident spectrum given in Eq. (9), if we consider only 
ionization losses and the range energy equation, R = aE r . At the present 
the above formulation will be terminated and the improvisions developed 
by the writer will be undertaken. 


The main improvement, by the writer is the introduction of an 
approximation for the proton range which represents the theoretical 
data, such as presented in Sternheimer 1 s article [l] , with an accuracy 
of ±4%, or better, for energies from around 5Mev to over 1,200 Mev. 
Also, the algebraic manipulation is essentially as elementary as that 
for the relationship, R = aE r . The new empirical formula for the range is 

R(E) = |- In (1 + 2bE r ) , (14) 

2b 

where a,b, and r are determined by fitting the range data of Ref. [ z] 
with the requirement to minimize the maximum relative error from 10 
to 1000 Mev. If, in Eq. (14), 2bE r << 1, then R = aE r . For example, 
in carbon, r = 1 . 78 , a = 2.3x10 3 , and b = 2 x 1 0 ^ , and we see that 
if E = 200 Mev, 2bE r = . 05 and using R =aE r one obtains an error to 
Eq. (14) of about 2. 5%. This good agreement is not obtained, however, 
with a larger Z number at such a large value of E . 

Figure 2 depicts an error analysis of the approximating function 
of Eq. (14) compared to data presented in Ref. 2 for two different co- 
efficients of r . In general practice it appears that for Z< 20, a value 
of r = 1. 78 is adequate, and for Z > 20, r— 1. 75 should be used. 
However, in the case of mixed materials of medium and low Z , it 
seems that a compromise may be made and that for a given calculation, 
one choice of r adhered to, perhaps 1. 78. Table I provides a summary 
of different values of a and b for different materials with r of 1 . 75 
and 1 . 78. It should be noted that a value of r = 1.8 is also given for 
tissue; this will be discussed in the development of the methods used 
by the writer for dose calculations. Figure 3 is a comparison of 
the error in the range for aluminum when using Eq. (14) to the error in 
range when using R = aE r . Because of the possible desire to employ 
other materials than those shown in Table I, the following relationships 
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FIGURE 2. THE RELATIVE ERROR IN APPROXIMATING PROTON RANGE USING EQ. (14) 





TABLE I 


Coefficients for the Range Equation 



r 

= 1. 75 


r = 

1. 78 

Material 

a 

b 


a 

b 

Carbon 

2 . 58x1 0~ 3 

1. 2x1 0" 6 


2 . 33x1 0 

2 . 0x1 0 6 

Aluminum 

3. 1 0x1 0 3 

1 . 9x10" 6 


2. 77x1 0 " 3 

2. 5x1 0" 6 

Iron 

3. 70x1 0 " 3 

2 . 6x1 0 6 


3 . 26x10 3 

3. 0x1 0" 6 

Copper 

3 . 85x1 0 3 

2. 7x1 0 " 6 


3 . 40x1 0 3 

3 . 25x1 0 

Silver 

Tungsten 

Polyethelene 

4. 55x1 0 3 

5. 50x1 0 3 
2 . 1 5x1 0 3 

3. 7x1 0" 6 
4. 2x1 0 ~ 6 
1 . lxlO' 6 


1 . 95x1 0 " 3 

1 . 7x1 0 ~ 6 

Tissue* 

2. 32x1 O -3 

1 . 2x1 0 " 6 


2. 1 lxlO' 3 

2 . 0x1 0 6 

Water 

2 . 32x1 0 3 

1 . 2x1 0 


2. 1 0x1 0 " 3 

2 . 0x1 0 6 

Air 

2 . 68x1 0 3 

1 . 4x1 0 -6 


2. 41xl0 -3 

2.1x10 6 

SiO z 

2. 87x1 0 3 

l . 7x1 0 ” 6 


2. 58x1 0 3 

2 . 5x1 0 6 

Glass 

3. 1 7x1 0 " 3 

2. lxlO' 6 


2. 83xl0" 3 

2 . 8x1 0 6 


_ - 3 

For stopping power in tissue: r Q = 1,80, a Q = 1.943x10 , b Q = 2. 273x10 . 
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were established by the writer and may be used for obtaining the co- 
efficients a and b when detailed curve fits are not warranted. 


r = 1.75 



1.6 x10 3 + 2. 89 X 10 - 4 YT 

5 . 1 6 x 1 0 " 7 V~Z 


(15) 


r = 1 . 78 


fa = 1. 53 x 10" 3 + 2. 33 x 10" ■ 4 Ya 

|b = 8x10 " 7 + 5x10 " 7 VT 


( 16 ) 


where A is the mass number and Z is the atomic number. 

The present calculations will be primarily limited to Z< Z0 ; 
hence, for Z < Z0 , we can write the range equation as 


R(E) 


1,53 x 10 3 + 2, 33 x 10 ~W]£ 
1.6 X 10"6 + 10 -6 iT 


In 1 4(1 . 6x 10" 6 VZ) E 1, 78 

(17) 


Even though the above equation is best for Z < 20, it is found that for 
practical shielding calculations the value of Z can be extended to about 
30 using r = 1 . 78, and the results are quite dependable if the shield 
thickness is greater than about two gm/cm 2 . The reason for this is 
that in copper (Z = 29), for example, the range of a 10-Mev proton is 
only about 0. 025 cm and a 1 0% error in range (. 0025 cm) at this ene rgy 
is trivial, if compared to the total shield thickness. It should be clear 
however that at higher energies the range error should be much smaller 
to maintain the above type of accuracy. 


Reverting to the original problem of this section, we develop 
the following relationships using Eq. (14) for the proton range. From 
Eq. (7), 


a 

2b 


In (l + 2bE r ) 


X + 2b ln 


(l + 2bE* r ) 


10 


or 


1 + 2bE r 
l + 2bE 


2bX 

a 


Solving for E, we obtain 


E = g(E*) 


*r V r 
(A+ BE ) 


where 


B 



and A = — (B - 1 ) . 
2b 


( 18 ) 


(19) 


From Eq. (19) it follows that 




and 


❖ 


= 0 if E < A 


Vt 


( 20 ) 


It is worth noting that if < < 1 , then A ^ X / a and B = 1 . 

(See Eq. 10.) For example, with carbon, 2bX/a = 1.717x10 3 X 
and for X^ 10 gm km 2 , the above approximation is quite valid. The 
foregoing analysis demonstrates why the simple range formula (R=aE r ) 
gives good results when X is not too large (X— 20 gm/cm 2 and 

‘l 

E< 2 50 Mev). Next, the differentiation of g(E ' ) gives 


dg(E v ) B E 

a = 

dE 


* r - 1 


r - 1 


(A+BE ) 


( 21 ) 


Substituting the above into Eq. (4), we obtain 


V E *> 


tg (eIE*)) b e* 1-1 

r - 1 / r 

(A + BE ) 


See discussion following Eq. (14). 


(2 2 ) 


1 1 


l 



There are two choices of the incident differential energy spectrum 
in vogue at present; the first is that given by Eq. (9) or else a family 
of N such curves given by 


<ME) 


Hi E 


-qi 


Ej < E 


^ E i + 1 


(23) 


where i = 1,2, 3, . . . , N. The second choice of representation is given 
by the integral rigidity spectrum 

- p / p ° 

N(>p) = N 0 e (protons/cm^); P > Pi (24) 


where p and p Q are in rigidity units of MV (millionvolts) . From 
Eq. (24) the differential rigidity spectrum becomes 


No " p /po 

<t> (p) dp = -dN(> p) = — e dp , p > pj . 


(25) 


In order to represent the above momentum rigidity units in 
energy (Mev) units, it is sufficient to use the relativistic relationship 
between variables given by (pze) 2 =E 2 4*2Em 0 1 or 

p = Ve 2 +1876E (26) 

where (ze) = 1 electron charge for protons, m 0 = 938 (the rest mass 
of the proton in Mev units), p is in MV and E is in Mev. Next, using 
a change of variable technique, we obtain 

P(E) 

0(E) dE = - dN (> p(E)) = ^2- e P ° 

Po 


dp 

dE 


dE , p> Pl , (27) 


1 Note that A V = work/q ; in basic physics, the potential difference 
is thus defined and, hence, Eq. (26) is dimensionally valid. 

2 m Q = 938.23 to 5 significant figures. 
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where 


dp = 

dE 


E + m. 


VE 7 T-"2m ”E" 


Substituting dp/dE intoEq. (27), we obtain 


1 /TT \ JTT - 
^ — 


N n (E + 938) 

P Q VE7TTT75TT 


:xt 


VE2 + 1876 E 


dE, E> Ej 

(28) 


where 


Ei = Vpi 2 + ( q 38) 2 - 938 , 


and <£ (E) has the units of protons/cm 2 Mev, The validity of the above 
transformation follows from elementary probability theory of distribution 
functions or else elementary calculus depending on the readers academic 
orientation. 


Referring to Eq. (22) it is of interest to obtain the proton differen- 
tial energy spectrum at a depth X using the incident spectrums of 
Eqs. (23) and (24). Using the incident spectrum of Eq. (23) we obtain: 


# X (E*) 


Hi B E 


* r r + qj - 1 
(A + B E'" ) r 




Vr 

(29) 


where Eq. (20) must be satisfied; B = exp(2bX/a) and A = (B-l)/2b. 
Using the rigidity spectrum of Eq. (24) we obtain (from Eq. 28): 




N 0 (s + 938) B E 


*r-l 


exp (- VP 


, r - 


I Vs^ + 18 76 s' 


+ 18 76 s/ p Q ) 


(30) 


where 


s 


(A + BE 



E > 



and Ei = Y Pl 2 +879,844 - 938. 
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The use of s was simply to shorten the size of the expression in Eq. 

(30). Examples of typical differential energy spectra as a function of 
depth X in aluminum are given in Figs. 4 and 5 illustrating the evalua- 
tion of Eqs. (29) and (30). 

In the following work, throughout this report, the presentation 
of proton penetration results (dose, flux, etc.) applies equally well to 
either plane slabs with a normal incident flux or at the center of spherical 
shells with an incident isotropic flux. This is due to the straight-ahead 
nature of the proton slowing down in a media, and the fact that the units 
of differential flux are protons/cm 2 - sec -Me v -steradian. Thus, if both 
time and direction are integrated out, the units are simply protons/cm 2 - 
Mev and the flux has no directional dependence. However, other informa- 
tion is usually provided concerning the source of the spectral data. If 
we are told that the spectra is for an omnidirectional flux, then the 
results apply only to the center of a sphere because of the (47 r) factor 
implicitly contained in the total flux. Otherwise the results are applicable 
as the reader sees fit. 


It is of some interest to note in Figs. 4 and 5 that the spectrum’s 
maximum shifts to the right as the slab thickness increases. It is a 
simple matter to find the value of E ' at which the spectrum's maximum 
occurs by solving the following equation for E^, 


9 gxjgQ 
a e"' 


o . 


(31) 


For example using Eq. (29) for 0^(E‘) one obtains 


(r-t) 


t- 


exp (•- 


2bX 


Vr 


Vr 


max 


2bqi 


(r-l)X 

aqi 


if 2*«,. 

“ (32) 


Equation (32) is plotted in Fig. 6 for three values of 
down media is aluminum with r - 1. 78, and 2b/a = 


q^. The slowing 
1 . 8 x 10' 3 . 
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FIGURE 4. PROTON DIFFERENTIAL ENERGY SPECTRUM AT 


DIFFERENT DEPTHS IN ALUMINUM 
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FIGURE 6. 



III. MULTILAYER SHIELDS 


The. above discussion is equally well applied to stratified layers 
of different materials. Figure 7 illustrates the parameters involved. 

In order to see the nature of the derivation for multiple layers of different 
materials, two layers will be considered initially. Referring to Eq. (19) 
let us define 


A i 


1 

2b x 



B x - exp 


ZbiX x 

a i 


and consequently, 


E o r = A ! + B ! E j r (33) 

where a! and hi are the material coefficients of Eq. (14); X x refers 
to the thickness of the first layer with E Q and Ej denoting the energies 
respectively incident on the first layer and transmitted through the first 
layer. Now applying the relationship of Eq. (7) to the second layer, 
we obtain 


ln {1 + 2b 2 E i r > = X 2 + Ur ln (1 + 2b 2 e/) ; 

l , L >2 ^02 


Simplifying, 


e* = a 2 + b 2 e 2 


(34) 


where 


B 2 - exp 


2b2 X 2 


and Az = 21 67 ' 0 


Substituting E x of Eq. (34) into Eq. (33), we obtain 


E Q r = Ai + B t (B 2 E 2 r + A 2 ) = (Aj + B,A 2 ) + B!B 2 E 2 r . (35) 


Equation (35) expresses the energy at a depth of X 2 in the 
second layer in terms of the energy incident on the first layer. If this 
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FIGURE 7. MULTILAYER SHIELDS AND ASSOCIATED PARAMETERS 



is repeated for N layers one obtains: 


r * ❖ >\c Y 

E = A + B E , (36) 

❖ 

where E is the incident energy on the first layer and E is the energy 
at the end of the N^ 1 layer; and, 

❖ 

A = Ax + A 2 B X + A 3 B x B 2 + ... + Ajsj B x B 2 B 3 . . . B]\j_| , 

(37) 

B" = BiB 2 B 3 . . . B n 


where B^ = exp | > Aj = (B^ - \ ) /Zb{ , and i = 1,2,..., N. 

This fairly simple representation of the energy as a function of depth 
and layer thicknesses of different materials is brought about by the 
fact that r is assumed to be constant for all materials considered. 

In shield optimization techniques, such a representation should be 
promising. Since Eq. (36) has the same form as Eq. (19), it follows 
that the coefficients A, B may be replaced by A' 1 ', B' whenever multi- 
layer shields are considered. Thus, all results obtained in the preceding 
or subsequent sections can be extended to multiple layers by using 
A*, B* for A,B. In the special case where 2bX/a << 1, i.e., (R^aE r ), 
then for the ith layer Bj = 1, A^ = Xj/ai and for N layers 


N 



It may be desirable at times to obtain an estimate of the proton 
transmission for two or more stratified layers of different materials 
but using only one material for the attenuation. Thus, it is necessary 
to find the equivalent thickness of the other layers in terms of the base 
material. 

This can be readily done in the following manner. First, we 
assume that for the base material, the simple power law holds for the 
range equation (R = aE r ). Then the thickness of the various layers in 
terms of the base element becomes simply, 
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where A denotes a mass or material number; a denotes the coefficient 
of the range equation 1 for the various materials; and the subscript \ 3 
denotes the base material. For example, if a reliable but simple 
estimate was desired for the proton dose rate behind 5 gm/cm 2 of 
Aluminum and 5 gms/cm 2 of tissue, the following information may be 
read off a plot of dose rates versus depth in aluminum by reading the 
dose rate at the depth: 


X 


A1 


X A1 + a" 


A1 


(tissue) 


X . 

tis sue 


or 



5 gm/cm 2 + 


2. 77 x 10 “ 
Z.llxlO -3 


(5 § m / cm 2 ) = 11.56 gm/cm 2 

(40) 


That the above technique is valid can be shown readily by inserting 
Xp into Eq. (38), thus obtaining 





The last equation denotes the equivalent energy transformation 
if we had used the material coefficients. Thus, the substitution of 
Eq. (39) gives the same results as the simple power law approximation 
of the range equations provided the same power r is used for all the 


1 The values of a are taken from Table I. 


materials. The utilization of the above simple relationships is of 
engineering importance for estimating depth dose due to primary protons 
when only a simple curve is provided for the so-called n skin dose. " 

See Section V on "Proton Dose Rate Calculations. " 


IV. NONELASTIC PROTON COLLISIONS AND SECONDARIES 

It was pointed out in the introduction to this paper that elastic 
scattering off a nucleus by high energy protons (>20 Mev) is highly 
forward with trivial reduction in energy. This assumption is not as 
valid for proton collisions in hydrogen but this problem will not be 
treated here. However, it is worth mentioning that the so-called range 
straggling associated with energetic protons is an effect mainly due 
to elastic collisions with electrons. However, this type of error is 
usually quite small and can be represented approximately for protons by 


o s 0. 015 R 
K 


(42) 


where a ^ is the standard deviation of a Gaussian distribution depicting 
the statistical fluctuation of the range about a mean range R (p. 662, 

[ 3 ]). This can be interpreted as meaning that 95% of monoener getic 
protons should have a measured range within about ± 3% of the theoretical 
range calculated from ionization losses only. This is not a bad error for 
shielding calculations since the proton energy spectrum always contains 
uncertainties of a much greater order of magnitude. This error is also 
in keeping with the use of the approximation for the range introduced by 
the writer (Eq. 14). Examination of the error curves in Fig. 2 shows 
that for energies from less than 1 0 Mev to over 1000 Mev the coefficients 
(a,b,r) can be chosen to maintain a maximum variation of less than 4% 
from an accurate theoretical calculation. 

In the treatment of nonelastic cross sections the writer has 
represented the cross section as a function of energy and mass number 
using an empirical expression which is amenable to obtaining closed 
form solutions in the mathematical operations which are necessary to 
obtain transmitted flux and dose rates. The greatest constraint in ob- 
taining an accurate expression for cross sections is the lack of adequate 
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experimental nonelastic cross sections in the range of 5 Mev to 50 Mev 
for protons. There are a few values at widely separated energies. 
However, the low-energy cross section seems to resemble that of 
neutrons to some extent and for energies from 5Mev to 18 Mev the 
nonelastic cross section of neutrons taken from Troubetzkoy [4 ] were 
used for the protons with a Coulomb correction in energy. Then the 
low-energy cross sections were blended into the proton nonelastic cross 
section at higher energies. For proton energies in the range of 200 to 
2000 Mev, the nonelastic cross section is fairly well represented by the 
theoretical formula of Fernbach, Serber, and Taylor, 


a 

ne 


7i R i 


l - (1 + 2KR)e 
2K 2 R 2 


■2KR 




cm 


(43) 


where R = r Q A 3 is the radius of the target nucleus with mass number 
A , and K" 1 is the mean free path in nuclear matter. A simpler ex- 
pression, determined by the author, which agrees well with experimental 
values, as well as Eq. (43), is given by 


a 


ne 


0. 38 



73 

[ barns ] 


(44) 


The reason for choosing the ratio (A/27) in Eq. (44) is that the non- 
elastic cross sections for aluminum (A = 27) will be the basis for the 
empirical formulas which are developed below. The requirements for 
such a formula are that the values of the cross section should be zero 
at zero energy, have a maximum between 5 and 25 Mev, and be approx- 
imately a constant (asymptote) as the energy exceeds say 200 Mev. 
Equation (45) satisfies these requirements with some degree of success, 
in addition.to being tailored for further mathematical operations: 


0. 38 
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ne 


(E) 




0-73 
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2r 


+ dE 


2r r 

E + fE + g 


(45) 


where d, f , g are constants to be determined by curve fitting techniques 
and r (= 1 . 78) is the same power as used in the range equation, (14). 
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Examination of Eq. (45) shows that cr ne (0) = 0. Also, as E becomes 
large, a ne (E) -*■ . 38 (A/27) -73 ; this is readily seen by dividing 

numerator and denominator of Eq. (45) by E^ r and letting E increase 
without bounds. In order to require that Eq. (45) has a proper maximum 
for some positive value of E = E M , the derivative of Eq. (45) is equated 
to zero and Eq. (46) is found; 


r _ eg + Vc‘ 2 g 2 - cdfg + d 2 g' 

M d - cf 


d > cf 


(46) 


73 

where c = . 38 (A/Z7)’ and the quantity under the radical is >0. 

It should be pointed out that Eq. (45) has a minimum for a negative 
value of E . In the variable E r , Eq. (45) is a n serpentine” (as en- 
countered in analytic geometry) which has been translated. In general, 
the type of curve which is represented by Eq. (45) is shown in Fig. 8. 



Figure 8. Nonelastic Cross Section Formula 

In order to determine the coefficients (d,f,g) of Eq. (45) it is 
sufficient to solve simultaneously the system of equations given by 

I) Xm d - cX M f - 2cX M g - dg = 0 
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II) 


(47) 


III) X Q d - a o X Q f - a o g = X Q 2 (a 0 .c) 


1 .78 


0 .73 


where X = E*”” and c = 0. 38 (A/27) ’ * . Equation I is derived 

from Eq. (46) and Eqs. II and III are from Eq. (45). When A=27 
(aluminum), the following coefficients were obtained: c = .38; d = 88; 

f = _85; g = 12, 000. The abo vs choice was determined by setting 
Em “ ^ 6 Mev with a M ” 1 barn and then requiring cr Q = .15 barn at 
5 Mev. The estimated values of o- ne (E < 20) for protons wer e 


found by using the neutron data of Troubetzkoy [ 4] 
relationship 


and applying the 


1 

CT p (E) = ff [E-B(z.A)] ; B(z,A) Sl.15 z/A 3 (48) 


where B(z, A) is the Coulomb potential barrier, E the incident proton 
energy, and a n is the total neutron nonelastic cross section taken from 
Ref. [4] for aluminum. Since for high energies (E > 200 Mev) the total 
nonelastic cross section is given by the simple formula of Eq. (44), the 
same sort of expression would be desirable for all energies. Hence the 
expression below is an attempt to derive such an expression: 


I 

ne 


(E, A) 


a (E) x 
ne 



CE 2r + dE r 
£2r + ffir + g 


(cm 2 /gm), 


(49) 


where 

N q = Avogadro’s number 
C = 8.479 X 10 -3 (27/A)° -27 , 

r = 1 . 78 , 

d = 1.9547 (A/27)- 221 _ - 27 = 1 . 9 5 4 7 (A/27) °‘ 12338 

f = -84. 9 (A/27)- 221 r = 84. 9 (A/27)- 39338 , and 

g = 1 1 , 996. 3 (A/27)- 442r = 12, 000 (A/27)- 78676 
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The writer makes no claim to. success in finding an adequate 
representation. However, the accuracy of Eq. (49) is probably sufficient 
for low-Z materials (Z < 30). Figure 9 shows the results of using the 
above fit for several different A numbers. Figure 10 is for a ne (E) 
in barns. The reversal of relative values between Fig. 9 and Fig. 10 
should be observed. 


In order to utilize the above equation in a computation, the 
following analysis is undertaken. If a proton energy flux 0 O (E) [protons/ 
cm 2 -Mev] is incident on a slab of material, A, and if nonelastic 
collisions are considered, the energy spectrum of the primary protons 
which get to a depth x(gm/cm 2 ) should be given by 



0 x( E *) e 



(x) dx 


(50) 


where 0 x (E ) is given by expressions such as Eqs. (29) and (30), and 
Z (x) is given by Eq. (49) with the energy E depending onxac- 
c<5hding to Eq. (19) or, simply, 


E 


(Qe 


l/x 


zb> 


Vi 


(51) 


where 


Q = 




and y 


2b 

a 


The problem hinges on the ability to integrate the exponential of Eq. (50). 
Hence the expression, 



(c - 2bd) - 4b (c - bd) Qe yx + 4c b z Q 2 e 2yX 

1 / XT 2 7/ V* 

1 - 2bf + b 2 g) - 4b ( 1 - bf) Qe + 4b 2 Q 2 e 


(52) 


is obtained by substituting the energy transformation of Eq. (51) into 
Eq. (49). Now, make the following substitutions: 
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FIGURE 9. PROTON NONELASTIC CROSS SECTIONS IN (CM 2 /GM) UNITS 





(barns) 



FIGURE 10. PROTON NONELASTIC CROSS SECTIONS IN BARNS/ATOM 
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ZbQe X = t and dx = (d is differential here) ; 

v t 

co — c - 2bd ; 0 = c - bd (d is a constant here) ; 
a - l - 2bf + 4b 2 g ; and j3 = 1 - bf . 


Then the integration of Eq. (52) can be represented by 




ibQe 


V x 


CO • 


20 1 + ct 2 


a - 2/3 1 + t2- 


dt 

7 ” 


(53) 


letting T -a - 2j3t + t 2 , we can write 

/ 

c_ I tdt 

z J 

2bQ 

Equation (54) is now a set of standard forms in most handbooks, 
one observation should be made about the expression 



(54) 


However , 


q = 4 a - 4/3 2 = 4b 2 (4g - f 2 ) , (55) 

where the q , as defined above, occurs in the handbook solutions to 
Eq. (54). Since the factor Vqp occurs in the handbook formulas, it 
should be noted that 


4g > f 2 (56) 

must be satisfied in order to obtain a valid solution. This inequality 
is readily satisfied as seen by examination of Eq. (49). Also, this 
inequality determines the nature of the solution since a solution is 
obtainable for Y-q but with a different mathematical form. Thus we 
obtain for the solution of Eq. (54), omitting details of simplification: 
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DX 4- log 


t 




* (x) dx = 


G 2 + (A +j + BE*V 

G 2 + (1 + E* r ) 2 


a 

4b 


(C-D) 


where 



D = 


C - 2bd 
l - 2bf + 4b 2 g 


G 



and 

F '= [f (D - C) - 4bg D ] . 


The definitions of A and B are given in Eq. (19). The value of F 
is approximately 3. 35 x 10” 4 for aluminum. Hence there appears to 
be some justification for dropping the arctan functions. That this 
assumption is valid will be shown below. In any case we can write 


r x 

-Jo S E *( x )dx 


, f * r 2 
G 2 + (j- + E r ) 

G 2 +(A +4+ BE" 1 ) 2 


4b <C - D) 


(58) 


x exp DX - F ^tan -1 


I A +— + BE 
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-| - tan | 


| +E *n 



of simply if F times the difference in the arctan functions is sufficiently 
small; 
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, ■ f * r 2 

G 2 + ( — + E r ) 

G 2 + (A + - + BE ) 2 
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^<C - D, 


e 


(59) 


Figure 11 is presented in order to illustrate the integral over penetra- 
tion depth of the differential nonelastic cross section of the proton slowing 
down from ionization losses and also the effect on the cross sections 
of the arctan terms in Eq. (57). The curves in Fig. 11 are simply the 
integral of Eq. (57) divided by the depth x of penetration in aluminum. 
These curves depict the energy dependence of the nonelastic cross 
sections and exemplifies dramatically the effects of the ionization losses. 
The dashed curves show the results obtained when the arctan terms 
are omitted. It should be noted that the cross sections are plotted as 
a function of the protons energy E‘ at the depth X in aluminum. 

Perhaps the next curve, Fig. 12, presents the most significant results 
of the foregoing analysis. Here Eq. (58) is plotted for the same thick- 
nesses as shown in Fig. 11. Rather interesting is the fact that the 
exponential function (Eq. 58 or 59) is very nearly a constant over the 
total energy range for a given penetration thickness. This is a very 
important demonstration since the examination of Fig. 11 might lead 
one to infer that the use of Eq. (57) or equivalent is necessary in order 
to obtain a reliable flux calculation for penetration thicknesses less 
than say 10 gm/cm 2 . Actually, the relative errors associated with 
using a constant cross section increase with depth so that at a depth of 
50 grams/cm 2 the error in the attenuation is of greater significance 
than at 2 gm/cm 2 . The major problem associated with choosing a 
constant cross section is the cross section dependence on a given proton 
energy spectrum. However, it is felt that an approximation such as 
given in Eq. (60) below is justified in view of the many uncertainties 
that exist in any energy spectrum or the total proton intensity. 



0 . 01 


(-! 

0.27 

cm 2 1 

\ a 
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QTQ 
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( 60 ) 


In Eq. (60) the constant, . 01, may be reduced to . 009 or less when 
X> 20 gm/cm 2 . However, the nonelastic cross section, as given by 
Eq. (60), should be restricted for application to the attenuation of the 
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FIGURE 11. AVERAGE PROTON NONELASTIC CROSS SECTION (CM 2 /GM) 
AT DIFFERENT DEPTHS IN ALUMINUM 
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FIGURE 12. ATTENUATION FACTOR OF PROT 




flux only. In order to estimate the production of secondary particles, 
Eq. (49) should be used. Thus, the intensity of nonelastic collisions 
per Mev at energy E' ,<: produced at the depth x is well approximated 

by the relationship: 


S X (E*) 



(E*,A) (j> x (E*) e 


. 009 



( 61 ) 


where the cross section E„(E*,A) is given by Eq. (49). The energy 
flux <£ X (E*) is given by either Eq. (29) or (30). From Eq. (6l), it 
follows that the total number of nonelastic collisions per gm at depth 
x is given by 


S 


x 



, _ . '!% 
S X (E'-') dE 


collisions 

gm 




(62) 


However, the results of Eq. (62) are not of importance since the quanti- 
ties desired are actually the number and energy of secondary neutrons 
or protons produced at depth x . A further discussion of secondaries 
will be treated in the next chapter on dose rates. In order to illustrate 
the nature of the curves depicted by Eq. (6l), an incident spectrum of 
the form 

- p /p 0 

N (> p) = N 0 e , p> Pi 


(Eqs. 24 and 30) was chosen. The results of this calculation are shown 
in Fig. 13. The curve in Fig. 14 illustrates the difference in shape 
between the flux curve and the nonelastic collision density curve. If 
a constant cross section was used for ^ ne (E ‘ , A) the curves in Fig. 14 
would be exactly parallel. 


V. PROTON DOSE RATE CALCULATIONS 


The next step in this development is to derive expressions for 
the primary proton tissue dose or dose rate as a function of shield 
thickness and/or depth dose in tissue. This is obtained as follows. 
For the general case after penetrating a depth x in a shield the dose 
rate is simply given by 
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E (mev) 


FIGURE 14. COMPARISON OF NONELASTIC PROTON COLLISIONS /GM -MEV TO PROTON 
DIFFERENTIAL ENERGY SPECTRUM IN 10 GM/CM 2 OF ALUMINUM 
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(63) 
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' S Z e( x ) dx 

<t> X (E*) S(E*) dE* 


where the energy E'" is taken at the penetration depth x . The S(E*) 
is the stopping power in tissue and is given by Eq. (1). See Eqs. (57) 
and (29) or (30) for the first two expressions in the integrand. The F 
is a flux-to-dose conversion factor depending on units of flux. The 
stopping power formula for tissue can be made compatible with the 
approximating range equation, (14), in the following manner. Using 
the definition of Eq. (2), we see that 


or 
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dE 
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dE 
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log e (1 +2b 0 E r °) 
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l o r o 


E l- r o + 2ba_ 
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(64) 


where a 0 ,b 0 ,r 0 are corresponding range coefficients for tissue (Fig. 15). 
Using the approximations suggested in the previous chapter for the non- 
elastic cross sections we can write the proton dose rate after trans- 
mitting several layers including tissue in the last layer in the following 
way: 
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(65) 


where 


F = 1.6x10 


-8 rads -cm 2 


, if HE*) has units of - p f°. ‘° ns , 

protons r cm 2 -Mev 
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or 


F = 5.76 x 10 


rads -cm 2 -sec 
hr -proton 


❖ 

if <j) (E ) has unit of 


protons 
cm 2 -sec -Mev 


The flux <j> (E *) is given by either Eq. (29) or (30) with the constants 
A*, B ^ defined for multiple layers as shown in Eq. (37). Also it 
should be noted that the r power used in Eq. (37) is constant for all 
layers; however, the r Q power used in the stopping power may be 
different. In fact in all calculations presented in this paper for dose 
the r is chosen to be 1. 78 for the shielding materials, but r Q is 
1.80 for the stopping power in tissue. This flexibility permits a small 
increase in accuracy with little loss in computational speed when numeri 
cal integration methods are employed. It should be pointed out that if 
Eq. (29) is used for the energy flux, then fqr each energy sector of the 
spectrum confined between two energies and E^ , another integral 

analogous to Eq. (65), is required, but the integration limits change 
with the H^,qi for each sector. However, this is conveniently carried 
out in a numerical integration process by using the coefficients H^qi 
which are necessary to satisfy the limits of Eq. (29) at the energy E*. 
Very often it is useful to examine the integrand as a function of E*. In 
this manner, a feeling is obtained for the important energy regions in 
terms of dose. Also, the slope of this curve should indicate the width 
of energy intervals necessary for an accurate numerical calculation. 
Thus, the differential proton dose is calculated as follows: 


dD x 

dE 


-(EjXi + EoX, + . . .) 

Fe 0(E*)S(E*) 


rads 

Mev 


( 66 ) 


Examples of Eq. (66) are shown in Fig. 16 and 17. The proton 
dose as a function of depth, calculated from 'Eq. (65), is shown in Figs. 
18 and 19. In the latter two figures, there is a curve labeled "Total 
Estimated Dose"; this refers to an approximation of dose in rads which 
corrects for the secondary protons and neutrons generated by nonelastic 
collisions. The correction is based on the observation that for a shield 
of low-Z materials the number of secondary protons and neutrons per 
nonelastic collision at energies below around 200 Mev is less than one 
cascade particle of each kind (protons and neutrons). With the above 
observation and other considerations, it became plausible to conjecture 
that if the nonelastic attenuation factor exp (-^ ne X) is omitted in the 
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TISSUE DOSE ( rad s/hr -mev) 



FIGURE 17. DIFFERENTIAL PROTON DOSE IN ALUMINUM 
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dose calculation, then a correction is made for the secondary particles. 
The foregoing is the correction made in this paper for the calculations 
of the "total Estimated Dose. ,! Thus, 


Total Est. Dose = Primary Proton Dose x exp (Z\Xi + E 2 X 2 + . . . ) 

(67) 

Of course such an approximation is only valid within certain fixed limits 
of shield thickness, Z number of target, and energy of colliding protons. 
However, to lend validity to the above assumption, Fig. 20 is presented. 
The secondary data in Fig. 20 were generated by C. W. Hill of Lockheed 
[2 ]. The interesting result is that the approximation of Eq. (67) is 
rather accurate for dose in rads for the thicknesses of aluminum shown. 
The approximation will probably become less dependable at greater 
thicknesses, but at these greater depths the total dose is substantially 
smaller and even a fairly large error in estimating secondary contri- 
butions may be unimportant from a practical point of view. Table II 
provides an analysis of Fig. 20. 


The numerical methods used in integrating Eq. (65) consisted 
simply of the following technique: 
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f (E) dE 
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£ 

V- 1 


£(kA - 



( 68 ) 


where 


A = 


b - a 
N 


Since the independent variable (E' 1 ') extends over such a large range, 
from 0 to over 1000 Mev, the summation of Eq. (68) is carried out in 
four separate sums so that the A can be increased in value as the 
energies increase. The number of terms in each sum is an input 
parameter which allows the user of the code to establish his own 
accuracy. 
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FIGURE 20. COMPARISON OF APPROXIMATION METHOD FOR TOTAL DOSE 
TO DETAILED CALCULATION OF TOTAL DOSE 







Some useful simplifications in Eq. (65) can be made if a power 
law input spectrum (Eq. 29) is used. Thus, the dose rate is represented 
by 
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* l / r 


E i+1 ~ A 

B'" 


Y* _ I 

HjB E 

r+a 

(A* + B*E* r ) * 

e'-A** 1 ^ 

B* 


*1 -r Q 


fa V1 - r 

i ~ 1 \ a o r o 


+ ^ dE*, 
a o r o I 

(69) 


where the A ,B are defined for multiple layers in Eq. (37). If we 
make the change of variables indicated by 


and 
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Now, making the substitutions, 


m = 


1+r -r ( 


n = 


r o^qj 


* r + 1 
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„* = 3Ljl± 
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( 72 ) 


we obtain, after some additional simplification. 
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(73) 


where the lower limit is set to zero if A*> E^ . This condition is 
met when the minimum proton energy has a range equal to or less 

than the minimum shield thickness. For example, if E^ = 30 Mev, then 
any aluminum thickness greater than 1.175 gm/cm 2 would cause the 
lower limit to be zero. Thus, one can always choose a thickness of 
shield so that the integral of Eq. (73) may be written as: 

D = <t> L J t m_1 (l-t) n_l dt + 2b 0 / t m *" 1 (l-t) n *“ l dt\ .(74) 

0 0 J 
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However, the integrals are now recognized as incomplete beta functions. 
Thus, one may write the dose as simply 


D = $ 


{< 


(// (m, n) + 2b 0 /3 a (m 


% , n*)^ 


(75) 


Now, if the assumption is made that the initial upper energy limit 

is sufficiently large, then the value of a approaches 1. Thus, 
the further simplification in terms of gamma functions is made: 
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r(m)r(n) 

r(m+n) 


+ 2b 


r(m*)r(n*) 


° r(m s >'+n*) 


-} 


q > 2 . (76) 


If the stopping power coefficient b 0 for tissue is set to zero we get 
simply: 


D = # 


r (m) r (n) 

r (m+n) 


q > 2 - r . 


(77) 


Equation (77) should be used when n*< 0, (q ^2). Finally, if all bi 
are set to zero for the range coefficients, then B*=l and A v = + 

X 2 /a 2 + . • . ; and if r Q = r we obtain the version of the simplest feasible 
model for proton dose rate calculations (See Eq. 13.), 
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FH exp (-giXi - E 2 X 2 - ...) 




a o r 
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l a i 
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(78) 


The above equations have the defect that the incident proton energy 
spectrum is represented by only one power function <f> o(E) = HE _< 1, E > E 0 
and in Eq. (78) the range is depicted by the simple relation R = aE r . 
However, the results obtained by using Eq. (76) are quite impressive 
as is demonstrated in Figs. 21, 22, and 23 by comparison to Alsmiller 
[5] and Hill [2] . In order to further validate the use of Eq. (76), the 
comparison of results using Eq. (76) and the same calculation using 
Eq. (65) with a numerical integration technique given by Eq. (68) is 
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FIGURE 21. COMPARISON OF ANALYTICAL CALCULATION OF 
DOSE TO NUMERICAL METHODS 
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TISSUE DOSE (rads) 



X (gm /cm*) 


FIGURE 22. COMPARISON OF ANALYTICAL CALCULATION 
OF DOSE TO NUMERICAL METHODS 
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made in Fig. 24. Figure 24 depicts the relative error between the 
gamma function computation (Dp) and the numerical method (D num ). 

For a shield thickness less than about 40 gm/cm 2 the Dp is about 
l/4% higher than the D num . This is probably due to the systematic 
error introduced in the numerical method which can be improved by 
taking more terms. (The numerical method used 55 steps between 0 
and 1000 Mev. ) However, as the thickness increases above 40 gm/cm 2 , 
the relative error increases more rapidly. This is due to two factors. 
The first is the fact that the gamma function represents an integration 
from an energy of zero to infinity. The second is the fact that the 
stopping power formula in tissue (Eq. 64) becomes a very poor approx- 
imation above 2000 Mev, resulting in a factor of six over estimation 
at 10, 000 Mev. One solution would be to refit the stopping power formula 
with more terms so that a more valid approximation is found for energies 
above 2000 Mev. However, it follows from consideratiqns of the proton 
differential spectrum that as E* exceeds 1000 Mev, the proton number 
approaches zero very rapidly. Hence, the results shown in Fig. 24 
imply that even for a thickness of 100 gm/cm 2 of aluminum the relative 
error due to the poor estimation of stopping power above 2000 Mev 
results in only a small over-estimate of dose depending fairly strongly 
on the power of the spectrum. For example, in Fig. 22, with q = -4.6 
the estimated error at 100 gm/cm 2 due to error in the stopping power 
above 1000 Mev leads to about a l/2% over estimate in dose, whereas 
when q = -3.95 this error is about 1%. The remaining differences 
can be explained by the integration to infinity for the gamma function, 
and the numerical integration error of about l/4% in Eq. (65). 

The foregoing analysis gave an indication of the numerical errors 
which occur in the proton dose rate calculation. However, a very im- 
portant source of error is often ignored in such discussions. What 
happens if the range energy data is systematically in error by some 
small amount? Analysis by the author has found that, in general, if 
the range is consistently low or high by, say, x%, then the dose is 
consistently low or high by about 2 x%. Thus a systematic error in range 
is reflected twice as great in the dose calculation. An example of the 
foregoing is presented in Fig. 25, and a detail error analysis of this 
figure and two other systematic errors are shown in Table III. 

The results of Table III, however, have one consoling aspect 
and this is that the errors are approximately uniformly displaced from 
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ANALYTICAL APPROXIMATION 






TISSUE DOSE (rads) 



FIGURE 25. EFFECT ON DOSE RATE AT CENTER OF SPHERE WHEN SYSTEMATIC 
ERROR OF 5% IS MADE IN PROTON RANGE 
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TABLE III 


Comparison of Relative Errors in Primary Proton Dose 
Where the Proton Range Has Been 
Systematically Under - or Over-Estimated 


gm/cm 2 

•s* 

Dose for 
Accurate 
Range 

Relative 
Error 
for +5%** 

Relative 
Error 
for - 5% 

1 

Relative 

Error 

for 4- 1 0% 

Relative 
Error 
for -10% 

Relative 
Error 
for + 1 5% 

Relative 
Error 
for -15% 

4. 5 

3. 33 lxl 0 3 

+.1010 

-.1104 

+. 2225 

-. 1988 

+ . 3425 

-. 2900 

9.0 

7. 38 7x1 0 2 

+.1010 

-.1105 

+. 2225 

-.1988 

+. 3427 

-. 2901 

18.0 

1 . 564x1 0 2 

+.1012 

-.1106 

+. 2228 

-. 199 1 

+. 3430 

-. 2905 

27. 0 

6. 0 72x1 0 1 

+. 1013 

-.1107 

+. 2231 

-.1993 

+ . 3435 

-. 2907 

54.0 

1 . 068x1 0 1 

+. 1018 

-.1112 

+. 2242 

-. 2C00 

+ . 3452 

-. 2922 

81 . 0 

3.437x10° 

+ . 1 026 

-.1119 

+. 2258 

-.2016 

+ . 3476 

-.2939 

108.0 

1.415x10° 

+. 1035 

-.1127 

+. 2276 

-.2031 

+ . 3505 

-. 2962 


Rads at center of sphere ; </> (E) - 6. 894xl0 15 E 3 ' 95 protons/cm 2 Mev. 

Relative error in dose when range is 5% too high but the stopping power in 
tissue is correct. 




TISSUE DOSE (rads) 



FIGURE 26. DOSE RATE AT CENTER OF SPHERE FROM A TYPICAL FLARE 
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the correct calculations. Thus one hopes that the deviations of the 
Bethe -Bloch range energy calculations vary randomly from the true 
proton range so that the error is periodic from say +5% to -5% and 
oscillating several times between 10 Mev and 1000 Mev, then the above 
type error in the dose calculation might well approach zero. 


As a final presentation in this section on Proton Dose Rates, 
a cursory analysis was made of solar proton events presented in Ref.[ 6] 
by W.R. Webber of the University of Minnesota. The result was a 
pseudo-average flare based on data from 1956 through 1962. This flare 
is represented by the integral rigidity spectrum 


N(>p) 


2. 5 x 10 9 



E 


protons 
cm2 -flare 


(79) 


where p is in units of MV (million volts), Rather interesting is the 
fact that even though the above flare is sort of an average, the probability 
is only 0. 025 that a more intense flare occurred during any one -week 
period between 1956 and 1962. Thus one could assume that the above 
flare represents a model for proton flares per week with a . 975 probability 
that the calculated doses will not be exceeded. The dose rates for the 
above spectrum at the center of a spherical enclosure of variable alum- 
inum thicknesses is presented in Fig. 26. 


VI. ISOTROPIC INCIDENT PROTON FLUX ON SLABS 


In the foregoing the calculations were for normal incident protons 
on plane slabs. Even though the calculations are denoted as dose rates 
at the center of spheres the computational geometry is for a monodirec- 
tional normal incident flux. If one wishes to extend the above calculations 
to oblique incidence, the only change necessary is that the slab layers 

be replaced by the thickness along the slant path of the oblique proton 
flux. Thus, the slant thickness is simply 


P { = X./cos 6 0 


(80) 
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where 6 Q is the angle between the slab normal and the angle of incidence 
of the proton flux. Thus it would be fairly simple, for any incident 
angular distribution, to calculate the dose rate at a given depth by cal- 
culating the proper solid angle weighting function and employing only 
the dose rate data for normal incidence. Thus, for the isotropic in- 
cident case one sees that if the normal incident dose, D(Xi), is calculated 
for the omnidirectional flux, then the primary proton dose rate for the 
isotropic incident flux on a slab is given by 


D 


Iso . 




X 

cos 6 


d cos 9 , 


or numerically the simple summation 


D 


Iso. 



(81) 


(82) 


where 


cos 0^ 




The above increment in cos 6 represents taking the i^h dose component 
at the midpoint of the ith solid angle. Also, the use of Eq. (82) would 
embody utilizing interpolation in a precalculated table of dose as a function 
of thickness. The above method was not utilized in the foregoing cal- 
culations even though it merits consideration as a computational method. 
The dose for an isotropic incident proton flux on a slab was calculated 
by employing numerical integration on the double integral below: 


D 


Iso. 



(E*)S(E*) sin 6 d0 dE* , (83) 


where the p= X/cos 0 and the calculations of B*,A* (Fig. 7) are 
calculated using = X^/cos 6 j for the slab penetration thicknesses. 

The results of the isotropic incident proton slab calculations are presented 
in Figs. 27 through 32. In Fig. 27, a comparison is made with Alsmiller’s 
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ISOTROPIC INCIDENT FLUX O N SLAB 

~4>(E ) « 2 .5 x id 1 ET 207 5 <E < 60 mev 

£(E)« 5.486 E' 3,98 60<E< 1000 mev 


mev-ster 


CALCULATED USING 2 « .01 cmVam 


CALCULATED BY ORNL 


X (gm/cm 2 ) OF AL 


FIGURE 27. SLAB PROTON DOSE RATE FROM ISOTROPIC INCIDENT FLUX 
COMPARED TO ALSMILLER [5] 






protons 




RADS/mev 



FIGURE 29. DIFFERENT PROTON DOSE AT DIFFERENT DEPTHS IN ALUMINUM 
FOR ISOTROPIC INCIDENT FLUX ON SLAB 


62 





TISSUE DOSE (RADS) 



FIGURE 30. PROTON DOSE AS FUNCTION OF SLAB THICKNESS FOR ISOTROPIC 
INCIDENT FLUX ON SLAB 
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TISSUE DOSE (rads) 



X (gm/cm ) OF AL 


FIGURE 31. COMPARISON OF DOSE FOR ISOTROPIC INCIDENT FLUX ON SLAB 
AND DOSE AT CENTER OF SPHERE 
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TISSUE (rads/steradlan) 



FIGURE 32. 


PROTON DOSE PER STERADIAN FOR ISOTROPIC INCIDENT FLUX 
ON S LAB OF ALUMINUM AT DIFFERENT DEPTHS 


TISSUE (rods/steradian) 



COS 9 


FIGURE 33. PROTON DOSE PER STERADIAN FOR ISOTROPIC INCIDENT FLUX 
ON SLAB OF ALUMINUM AT DIFFERENT DEPTHS 




work [5] . The agreement seems rather good. Figures 28 and 29 
are presented to parallel the treatments given in the previous chapters 
of this report for the same incident rigidity spectrum. Of course, the 
differential flux and dose/Mev are defined the same as previously ex- 
cepting that an integration was performed over direction for the iso- 
tropic slab case. Figure 30 presents the double integration of Eq. (83) 
for the same incident flux as shown in Figs. 28 and 29. Figure 31 is 
presented merely to emphasize the difference between two geometries 
and compare the primary dose of Fig. 18 with that of Fig. 30. It should 
be noted that at X = 0 gm/cm 2 depth, the isotropic dose on a slab is 
l/2 the dose at the center of a sphere. Thus, there is a systematic 
difference of a factor of two due to the difference of a 2 ir space and 
a 4 7r space in the inferred solid angle integration. Figures 32 and 33 
complete the isotropic incidence analysis. Here the dose/steradian is 
presented for several thicknesses of aluminum ranging from 1 gm/cm 2 
to 50 gm/cm 2 in Fig. 33. The possible applications of this type data 
are probably small since the straight-ahead model used in proton pene- 
tration calculations permits fairly simple numerical methods for even 
complex geometries. However, the curves do indicate the relative 
significance of the direction of incidence on the primary proton dose 
in finite slabs. The integration over solid angle of these curves gives 
the data points shown in Fig. 30. The curves in Figs. 32 and 33 indicate 
clearly that care should be exercised in an arbitrary assumption such 
as any proton with an angle of incidence greater than 45 to the normal 
can be ignored. For example, even for 10 gm/cm 2 , one would probably 
feel it necessary to consider angles as large as 70 (cos 8 - . 34) in 
order to obtain a reliable integration over angle. The implications of 
Figs. 32 and 3 3 on the treatment of complex geometries leads one to 
recognize that relatively large proton fluxes may enter a detector at 
rather oblique angles measured from the vehicle surface normal, 
particularly if the vehicle walls are relatively thin. A study of more 
complex geometry effects is given by the writer in the report "Flare 
Proton Doses Inside Lunar Structures” [ 7 ] . Three examples taken 
from this report are shown in Figs. 34 and 36. The incident proton 
flux in these calculations was assumed to be isotropic and the rigidity 
spectrum N(>p) = 4. 54 x 10 10 e used throughout this report 

was assumed. The dose calculations in Figs. 34, 35, and 36 are based 
on the "Total Estimated Dose,” which include s secondaries. 
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FIGURE 34. CENTER LINE DOSE VS. DISTANCE FROM BASE OF FLAT 
TOP CYLINDER 
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FIGURE 35. CENTER LINE DOSE VS. DISTANCE FROM BASE OF CYLINDER 
WITH SPHERICAL CAP 
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DOSE (RADS) 



FIGURE 36. AXIS DOSE VS. DISTANCE FROM GROUND OF HEMISPHERE AND SPHERE 
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VII. CONCLUSIONS 


One of the difficulties in an analysis of the type presented in 
the preceding pages is finding a point at which to stop. There are several 
areas which have been mentioned but not subjected to adequate study. 

For example, the examination of multiple layer shields and depth dose 
in tissue, as well as the tedious problem of treating secondary particle 
production fractions and their attenuation in a shield. However, it is 
hoped that usable short cuts and insight into various aspects of proton 
shielding have been made. Also, it is hoped that adequate information 
has been presented in the report for determining the validity of the 
results obtained and the scope of applications of the suggested approxi- 
mations . 

The computer codes used in the foregoing work are written in 
Fortran and are available to anyone interested. The proton dose rate 
code is capable of treating up to ten layers of different materials and 
takes directly either the coefficients of the integral rigidity spectrum 
(P in Million volts) or the power law representation of the differential 
spectrum (protons/cm 2 -Mev-sec). The power law spectrum can be 
broken into as many as ten segments or energy groups. The output is 
rather extensive depending on the users’ needs. Thus, differential 
flux, dose, and nonelastic collision density can be found as a function 
of the proton's energy inside the shield. If a choice is made to treat 
isotropic incidence of protons on a slab, the angular distribution of the 
transmitted proton dose is printed out for up to 20 increments in the 
angle Q measured from the slab normal. The greatest virtue of the 
above codes is the lack of large amounts of input data since the approxi- 
mations discussed in this report are utilized wherever there is no real 
compromise in accuracy. 
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